##########################################################################################

library(dplyr)
library(ggplot2)
library(data.table)
library(RColorBrewer)
library(optparse)
library(ggpubr)

##########################################################################################

option_list <- list(
    make_option(c("--gene"), type = "character") ,
    make_option(c("--info_file"), type = "character") ,
    make_option(c("--mut_rate_gene_file"), type = "character") ,
    make_option(c("--mut_rate_point_file"), type = "character") ,
    make_option(c("--images_path"), type = "character")
)

if(1!=1){
    
    gene <- "TP53"
    work_dir <- "~/20220915_gastric_multiple/dna_combinePublic/"
    info_file <- paste(work_dir,"/config/tumor_normal.class.list",sep="")
    mut_rate_gene_file <- paste(work_dir,"/images/mutRate/MutRate.IM.tsv",sep="")
    mut_rate_point_file <- paste(work_dir,"/images/mutRate/MutRate.RecurrentPoint.IM.tsv",sep="")
    images_path <- paste(work_dir,"/images/mutRatePlot",sep="")

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

gene <- opt$gene
info_file <- opt$info_file
mut_rate_gene_file <- opt$mut_rate_gene_file
mut_rate_point_file <- opt$mut_rate_point_file
images_path <- opt$images_path

dir.create(images_path , recursive = T)

###########################################################################################

col_im <- c(brewer.pal(9,"YlGnBu")[6:8])
names(col_im) <- c("IM(IGC)" , "IM(DGC)" , "IM(IGC_DGC)")
col_im <- col_im[1:2]
col <- col_im[1:2]

###########################################################################################

dat_mutRateGene <- data.frame(fread( mut_rate_gene_file ))
dat_mutRatePoint <- data.frame(fread( mut_rate_point_file ))
dat_info <- data.frame(fread( info_file ))

###########################################################################################
## 计算每种亚型IM的数量
dat_sampleNum <- dat_info %>%
group_by( Type ) %>%
summarize( SampleNum = length(unique(ID)) )

dat_sampleNum <- data.frame( dat_sampleNum )
dat_sampleNum$Type <- ifelse( dat_sampleNum$Type == "IM + IGC" , "IM(IGC)" ,dat_sampleNum$Type )
dat_sampleNum$Type <- ifelse( dat_sampleNum$Type == "IM + DGC" , "IM(DGC)" ,dat_sampleNum$Type )
dat_sampleNum$Type <- ifelse( dat_sampleNum$Type == "IM + IGC + DGC" , "IM(IGC_DGC)" ,dat_sampleNum$Type )

###########################################################################################

dat_mutRateGene <- subset(dat_mutRateGene , Hugo_Symbol==gene)
dat_mutRatePoint <- subset(dat_mutRatePoint , Hugo_Symbol==gene)
dat_mutRatePoint$Hugo_Symbol <- dat_mutRatePoint$vid
dat_mutRatePoint <- subset( dat_mutRatePoint , vid %in% unique(subset(dat_mutRatePoint , MutNum > 5)$vid) )
dat_mutRatePoint <- dat_mutRatePoint[,-3]

dat_plot <- rbind( dat_mutRateGene , dat_mutRatePoint )
dat_plot$Class <- factor( dat_plot$Class , levels = names(col_im) , order = T )
dat_plot$value_text <- paste0( round(dat_plot$MutRate , 2) * 100 , "%") 
dat_plot <- subset(dat_plot , Class != "IM(IGC_DGC)")

###########################################################################################

trans <- function(num){
    up <- floor(log10(num))
    down <- round(num*10^(-up),2)
    text <- paste("p == ",down," %*% 10","^",up)
    return(text)
}

###########################################################################################
## 计算P值
dat_plot$p.value = ""
dat_plot$p_text = ""
#dat_plot <- subset( dat_plot , Hugo_Symbol %in% unique(dat_plot[dat_plot$MutNum > 2,"Hugo_Symbol"]))

class_compare <- names(col_im)

result <- c()
for(geneN in unique(dat_plot$Hugo_Symbol)){

    print(geneN)

    for( i in 1:(length(class_compare)-1) ){
        class1 <- class_compare[i]
        for( j in (i+1):length(class_compare) ){
            class2 <- class_compare[j]
            tmp_1 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% class1 )
            tmp_2 <- subset( dat_plot , Hugo_Symbol == geneN & Class %in% class2 )

            if(nrow(tmp_1)==0){
                tmp_1 <- tmp_2
                tmp_1$SampleNum <- dat_sampleNum[dat_sampleNum$Type==class1,"SampleNum"]
                tmp_1$MutNum <- 0
                tmp_1$MutRate <- 0
                tmp_1$value_text <- ""
            }

            if(nrow(tmp_2)==0){
                tmp_2 <- tmp_1
                tmp_2$SampleNum <- dat_sampleNum[dat_sampleNum$Type==class2,"SampleNum"]
                tmp_2$Class <- class2
                tmp_2$MutNum <- 0
                tmp_2$MutRate <- 0
                tmp_2$value_text <- ""
            }

            tmp <- rbind( tmp_1 , tmp_2 )
            tmp_fisher <- matrix(c(tmp$MutNum , tmp$SampleNum - tmp$MutNum) , ncol = 2)
            p <- fisher.test(tmp_fisher)$p.value

            if( p < 0.001 ){
                p_text <- trans(p)
            }else{
                p_text <- paste0( "p == " , round(as.numeric(p) , 3) ) 
            }
           
            tmp$p.value <- p
            tmp$p_text <- p_text

            result <- rbind( result , tmp )
        }
    }
}

result$Hugo_Symbol <- factor( result$Hugo_Symbol , 
    levels = unique(result[order(result$MutRate , decreasing=T),"Hugo_Symbol"]) , order = T)
result$Class <- factor( result$Class , levels = names(col_im) , order = T )
out_name <- paste0( images_path , "/MutRate_" , gene , ".IM.tsv" )
write.table( result , out_name , row.names = F , sep = "\t" , quote = F )

###########################################################################################

dat_use <- unique( result[,1:6] )

plot <- ggplot( data = dat_use , aes( x = Class , y = MutRate , fill = Class ))+
geom_bar(position = "stack", stat = "identity") + 
theme_bw()+
labs(x="",y="Mutation Rate")+
facet_grid(.~Hugo_Symbol) +
theme(panel.grid = element_blank())+
scale_fill_manual(values=col_im) +
ylim(0,1)+
geom_text(aes(label=value_text) , position=position_stack(vjust = 0.5) , size=2.5 , color="black")+
geom_text(aes(label=p_text , y = 1 ,x = 1.5),parse = TRUE,size=4)+
theme(panel.background = element_blank(),#设置背影为白色#清除网格线
            legend.position ='none',
            legend.title = element_blank() ,
            panel.grid.major=element_line(colour=NA),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 7,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            strip.text.x = element_text(size = 7 , face = 'bold'),
            axis.ticks.x = element_blank(),
            axis.text.x = element_text(size = 8,color="black",face='bold') ,
            axis.line = element_line(size = 0.5))

out_name <- paste0( images_path , "/MutRate_" , gene , ".IM.pdf" )
ggsave( out_name , plot , width = (2 * length(unique(dat_plot$Hugo_Symbol))) , height = 4 )
